function [fishMat, snr, AeffWeightd] = computefishermatrix_ifo(Params, psdType, lambdaVec)
%
% COMPUTEFISHERMATRIX_IFO - compute the fisher information matrix from the 3.5PN 
% accurate restricted PN waveforms in the SPA. 
%
% usage: [fishMat, snr]=computefishermatrix_ifo(Params, psdType, lambdaVec)
%
% Params    : structure containing parameters totalMass and eta 
% psdType   : type of noise curve to use: LIGO/AdvLIGO/Virgo/WhiteNoise
% lambdaVec : a cell array of parameters with which the Fisher matrix should be computed. 
% 
% Anand Sengupta, 13.03.2008
%
% $Id: computefishermatrix_ifo.m 68 2010-01-20 21:58:45Z anand.sengupta $

 setconstants

 P = Params.P;
 Q = Params.Q;

 fLow = Params.fLower;
 flso = (1./sqrt(6))^3 / (pi*Params.totalMass*MSOLAR_TIME);

 if flso <= fLow 
 	fishMat = [];
	warning('Fisher matrix is empty. flso < fLower')
	return
 end

 df   = 2.5e-1;
 fVec = [fLow:df:flso];

 % generate noise PSD
 [Shf] = detnoisepsd(psdType,fVec);

 % compute SNR 
 Aeff		= sqrt(Params.P^2 + Params.Q^2);
 AeffWeightd	= Aeff * fVec.^(-7/6) .* sqrt(df./Shf);
 snr 		= 2 * Aeff * sqrt(dot(fVec.^(-7/3), 1./Shf) *df);

 % compute each element of the fisher matrix
 for i=1:length(lambdaVec)
	for j=1:length(lambdaVec)

	   % compute the derivatives of the phase and amplitude with respect
	   % to parameter lambda_i
	   [dPsieffBydLambda_i]	= PsieffDerivative(lambdaVec{i}, fVec, Params);
	    dPBydLambda_i		= PDerivatives(lambdaVec{i}, Params);
	    dQBydLambda_i		= QDerivatives(lambdaVec{i}, Params);

	   % compute the derivatives of the phase and amplitude with respect
	   % to parameter lambda_j
	   [dPsieffBydLambda_j]	= PsieffDerivative(lambdaVec{j}, fVec, Params);
	    dPBydLambda_j		= PDerivatives(lambdaVec{j}, Params);
	    dQBydLambda_j		= QDerivatives(lambdaVec{j}, Params);

	   % now compute the fisher matrix element by taking the noise-weighted
	   % inner products. See the paramEstNotes document for clarification.

	   dHidHj = (dPBydLambda_i*dPBydLambda_j + dQBydLambda_i*dQBydLambda_j).*fVec.^(-7/3);
	   dHidHj = dHidHj...
		 + (P^2 + Q^2)*dPsieffBydLambda_i.*dPsieffBydLambda_j.*fVec.^(-7/3);

	   fishMat(i,j) = 4 * dot(dHidHj, 1./Shf) * df;

	end
 end

 return;








